function image_compression_svd()
% IMAGE_COMPRESSION_SVD 基于SVD的图像压缩专项实验
% 
% 详细研究SVD在图像压缩中的应用和优化技术

fprintf('=== 基于SVD的图像压缩专项实验 ===\n\n');

%% 实验1：基础图像压缩
fprintf('实验1：基础图像压缩\n');
fprintf('目标：理解SVD图像压缩的基本原理\n\n');

% 创建测试图像
img_size = 128;
[X, Y] = meshgrid(1:img_size, 1:img_size);

% 合成图像：包含不同频率成分
test_image = 0.5 * sin(2*pi*X/img_size) .* cos(2*pi*Y/img_size) + ...
             0.3 * sin(4*pi*X/img_size) + ...
             0.2 * cos(6*pi*Y/img_size) + ...
             0.1 * randn(img_size, img_size);

% 归一化到[0,1]
test_image = (test_image - min(test_image(:))) / (max(test_image(:)) - min(test_image(:)));

fprintf('测试图像尺寸: %d×%d\n', size(test_image));
fprintf('像素值范围: [%.3f, %.3f]\n', min(test_image(:)), max(test_image(:)));

% SVD分解
[U1, S1, V1] = svd(test_image);
singular_values = diag(S1);

fprintf('奇异值统计:\n');
fprintf('最大奇异值: %.4f\n', singular_values(1));
fprintf('最小奇异值: %.4f\n', singular_values(end));
fprintf('条件数: %.2e\n', singular_values(1) / singular_values(end));

% 奇异值能量分布
energy_cumulative = cumsum(singular_values.^2) / sum(singular_values.^2);
fprintf('前10个奇异值的累积能量占比:\n');
for i = 1:10
    fprintf('前%2d个: %6.2f%%\n', i, 100 * energy_cumulative(i));
end

%% 实验2：压缩质量分析
fprintf('\n实验2：压缩质量分析\n');
fprintf('目标：分析不同压缩级别的质量指标\n\n');

compression_ranks = [1, 2, 5, 10, 20, 30, 50, 80];
fprintf('秩    压缩比    PSNR(dB)    SSIM      MSE       存储节省\n');
fprintf('--    ------    --------    ----      ---       --------\n');

% 存储压缩结果用于可视化
compressed_images = cell(length(compression_ranks), 1);

for i = 1:length(compression_ranks)
    r = compression_ranks(i);
    
    % SVD压缩重构
    img_compressed = U1(:, 1:r) * S1(1:r, 1:r) * V1(:, 1:r)';
    compressed_images{i} = img_compressed;
    
    % 压缩比计算
    original_elements = numel(test_image);
    compressed_elements = r * (size(test_image, 1) + size(test_image, 2) + 1);
    compression_ratio = compressed_elements / original_elements;
    storage_saving = (1 - compression_ratio) * 100;
    
    % 质量指标
    mse = mean((test_image(:) - img_compressed(:)).^2);
    psnr = 10 * log10(1 / mse);
    
    % SSIM计算（简化版本）
    ssim_val = compute_ssim(test_image, img_compressed);
    
    fprintf('%2d    %6.3f    %8.2f    %6.4f    %.2e    %7.1f%%\n', ...
            r, compression_ratio, psnr, ssim_val, mse, storage_saving);
end

% 可视化压缩结果
figure('Name', 'SVD图像压缩结果', 'Position', [100, 100, 1200, 800]);
subplot_positions = [1, 2, 3, 4, 5, 6, 7, 8];
for i = 1:min(8, length(compression_ranks))
    subplot(2, 4, subplot_positions(i));
    imagesc(compressed_images{i});
    colormap gray; axis equal; axis off;
    
    r = compression_ranks(i);
    mse = mean((test_image(:) - compressed_images{i}(:)).^2);
    psnr = 10 * log10(1 / mse);
    title(sprintf('秩 %d (PSNR: %.1f dB)', r, psnr));
end

%% 实验3：不同图像类型的压缩效果
fprintf('\n实验3：不同图像类型的压缩效果\n');
fprintf('目标：比较不同类型图像的SVD压缩效果\n\n');

% 创建不同类型的测试图像
image_types = {
    '平滑渐变', @(x,y) (x/img_size + y/img_size) / 2;
    '棋盘图案', @(x,y) mod(floor(x/16) + floor(y/16), 2);
    '随机噪声', @(x,y) rand(size(x));
    '环形图案', @(x,y) exp(-((x-img_size/2).^2 + (y-img_size/2).^2) / (img_size/4)^2);
    '条纹图案', @(x,y) sin(8*pi*x/img_size);
};

fprintf('图像类型      秩10压缩比    PSNR(dB)    有效秩估计\n');
fprintf('--------      ----------    --------    ----------\n');

for i = 1:length(image_types)
    type_name = image_types{i, 1};
    img_func = image_types{i, 2};
    
    % 生成图像
    test_img = img_func(X, Y);
    test_img = (test_img - min(test_img(:))) / (max(test_img(:)) - min(test_img(:)));
    
    % SVD分解
    [U_type, S_type, V_type] = svd(test_img);
    s_type = diag(S_type);
    
    % 秩10压缩
    r_test = 10;
    img_comp = U_type(:, 1:r_test) * S_type(1:r_test, 1:r_test) * V_type(:, 1:r_test)';
    
    % 质量指标
    compression_ratio = (r_test * (img_size * 2 + 1)) / (img_size^2);
    mse = mean((test_img(:) - img_comp(:)).^2);
    psnr = 10 * log10(1 / mse);
    
    % 有效秩估计（90%能量）
    energy_cum = cumsum(s_type.^2) / sum(s_type.^2);
    effective_rank = find(energy_cum >= 0.9, 1);
    
    fprintf('%-12s  %10.3f    %8.2f    %10d\n', ...
            type_name, compression_ratio, psnr, effective_rank);
end

%% 实验4：自适应秩选择
fprintf('\n实验4：自适应秩选择\n');
fprintf('目标：开发自动选择最优压缩秩的方法\n\n');

% 使用原始测试图像
[U4, S4, V4] = svd(test_image);
s4 = diag(S4);

% 方法1：能量阈值法
energy_thresholds = [0.90, 0.95, 0.99, 0.999];
fprintf('自适应秩选择方法比较:\n');
fprintf('方法                    选择的秩    PSNR(dB)    压缩比\n');
fprintf('----                    --------    --------    ------\n');

for threshold = energy_thresholds
    energy_cum = cumsum(s4.^2) / sum(s4.^2);
    optimal_rank = find(energy_cum >= threshold, 1);
    
    % 重构图像
    img_adaptive = U4(:, 1:optimal_rank) * S4(1:optimal_rank, 1:optimal_rank) * V4(:, 1:optimal_rank)';
    
    % 质量评估
    mse_adaptive = mean((test_image(:) - img_adaptive(:)).^2);
    psnr_adaptive = 10 * log10(1 / mse_adaptive);
    compression_ratio_adaptive = (optimal_rank * (img_size * 2 + 1)) / (img_size^2);
    
    fprintf('能量阈值 %.1f%%         %8d    %8.2f    %6.3f\n', ...
            threshold*100, optimal_rank, psnr_adaptive, compression_ratio_adaptive);
end

% 方法2：奇异值跳跃检测
fprintf('\n奇异值跳跃检测:\n');
s_ratios = s4(1:end-1) ./ s4(2:end);
[max_ratio, jump_idx] = max(s_ratios);
fprintf('最大奇异值比值: %.2f (位置: %d)\n', max_ratio, jump_idx);

img_jump = U4(:, 1:jump_idx) * S4(1:jump_idx, 1:jump_idx) * V4(:, 1:jump_idx)';
mse_jump = mean((test_image(:) - img_jump(:)).^2);
psnr_jump = 10 * log10(1 / mse_jump);
compression_ratio_jump = (jump_idx * (img_size * 2 + 1)) / (img_size^2);

fprintf('跳跃检测法              %8d    %8.2f    %6.3f\n', ...
        jump_idx, psnr_jump, compression_ratio_jump);

%% 实验5：分块SVD压缩
fprintf('\n实验5：分块SVD压缩\n');
fprintf('目标：使用分块策略提高压缩效率\n\n');

% 分块参数
block_sizes = [16, 32, 64];
target_rank = 5;

fprintf('分块大小    块数量    平均PSNR    总压缩比    计算时间\n');
fprintf('--------    ------    --------    --------    --------\n');

for block_size = block_sizes
    tic;
    
    % 分块处理
    [m, n] = size(test_image);
    blocks_m = ceil(m / block_size);
    blocks_n = ceil(n / block_size);
    total_blocks = blocks_m * blocks_n;
    
    img_blocked = zeros(size(test_image));
    psnr_blocks = [];
    
    for i = 1:blocks_m
        for j = 1:blocks_n
            % 提取块
            row_start = (i-1) * block_size + 1;
            row_end = min(i * block_size, m);
            col_start = (j-1) * block_size + 1;
            col_end = min(j * block_size, n);
            
            block = test_image(row_start:row_end, col_start:col_end);
            
            % 块SVD压缩
            [U_block, S_block, V_block] = svd(block);
            r_block = min(target_rank, min(size(block)));
            block_compressed = U_block(:, 1:r_block) * S_block(1:r_block, 1:r_block) * V_block(:, 1:r_block)';
            
            % 存储压缩块
            img_blocked(row_start:row_end, col_start:col_end) = block_compressed;
            
            % 块PSNR
            block_mse = mean((block(:) - block_compressed(:)).^2);
            if block_mse > 0
                block_psnr = 10 * log10(1 / block_mse);
                psnr_blocks = [psnr_blocks; block_psnr];
            end
        end
    end
    
    compute_time = toc;
    
    % 整体质量评估
    overall_mse = mean((test_image(:) - img_blocked(:)).^2);
    overall_psnr = 10 * log10(1 / overall_mse);
    avg_block_psnr = mean(psnr_blocks);
    
    % 压缩比估算
    elements_per_block = block_size^2;
    compressed_elements_per_block = target_rank * (block_size * 2 + 1);
    block_compression_ratio = compressed_elements_per_block / elements_per_block;
    
    fprintf('%6d      %6d    %8.2f    %8.3f    %8.4f\n', ...
            block_size, total_blocks, avg_block_psnr, block_compression_ratio, compute_time);
end

%% 实验6：颜色图像压缩
fprintf('\n实验6：颜色图像压缩\n');
fprintf('目标：扩展SVD压缩到彩色图像\n\n');

% 创建简单的彩色图像（RGB三通道）
color_img = zeros(64, 64, 3);
[X_color, Y_color] = meshgrid(1:64, 1:64);

% 红色通道
color_img(:, :, 1) = sin(2*pi*X_color/64) .* cos(2*pi*Y_color/64);
% 绿色通道
color_img(:, :, 2) = cos(3*pi*X_color/64);
% 蓝色通道
color_img(:, :, 3) = sin(4*pi*Y_color/64);

% 归一化
for c = 1:3
    channel = color_img(:, :, c);
    color_img(:, :, c) = (channel - min(channel(:))) / (max(channel(:)) - min(channel(:)));
end

fprintf('彩色图像尺寸: %d×%d×%d\n', size(color_img));

% 分通道SVD压缩
ranks_color = [1, 3, 5, 10];
fprintf('压缩秩    红通道PSNR    绿通道PSNR    蓝通道PSNR    总体PSNR\n');
fprintf('------    ----------    ----------    ----------    --------\n');

for r = ranks_color
    color_compressed = zeros(size(color_img));
    channel_psnrs = zeros(3, 1);
    
    for c = 1:3
        % 单通道SVD
        [U_c, S_c, V_c] = svd(color_img(:, :, c));
        color_compressed(:, :, c) = U_c(:, 1:r) * S_c(1:r, 1:r) * V_c(:, 1:r)';
        
        % 通道PSNR
        channel_mse = mean((color_img(:, :, c) - color_compressed(:, :, c)).^2, 'all');
        channel_psnrs(c) = 10 * log10(1 / channel_mse);
    end
    
    % 总体PSNR
    overall_mse = mean((color_img(:) - color_compressed(:)).^2);
    overall_psnr = 10 * log10(1 / overall_mse);
    
    fprintf('%4d      %10.2f    %10.2f    %10.2f    %8.2f\n', ...
            r, channel_psnrs(1), channel_psnrs(2), channel_psnrs(3), overall_psnr);
end

%% 实验7：压缩算法比较
fprintf('\n实验7：压缩算法比较\n');
fprintf('目标：比较SVD与其他压缩方法\n\n');

% 实现简单的DCT压缩作为对比
function img_dct = dct_compress(img, keep_ratio)
    % 简化的DCT压缩
    img_dct_full = dct2(img);
    
    % 保留左上角的系数
    [m, n] = size(img_dct_full);
    keep_m = round(m * sqrt(keep_ratio));
    keep_n = round(n * sqrt(keep_ratio));
    
    img_dct_compressed = zeros(size(img_dct_full));
    img_dct_compressed(1:keep_m, 1:keep_n) = img_dct_full(1:keep_m, 1:keep_n);
    
    img_dct = idct2(img_dct_compressed);
end

% 比较不同压缩方法
compression_ratios_test = [0.1, 0.2, 0.3, 0.5];
fprintf('压缩比    SVD PSNR    DCT PSNR    SVD优势\n');
fprintf('------    --------    --------    -------\n');

for ratio = compression_ratios_test
    % SVD压缩
    target_elements = round(ratio * numel(test_image));
    % 估算需要的秩
    estimated_rank = round(target_elements / (img_size * 2 + 1));
    estimated_rank = max(1, min(estimated_rank, min(size(test_image))));
    
    img_svd = U1(:, 1:estimated_rank) * S1(1:estimated_rank, 1:estimated_rank) * V1(:, 1:estimated_rank)';
    mse_svd = mean((test_image(:) - img_svd(:)).^2);
    psnr_svd = 10 * log10(1 / mse_svd);
    
    % DCT压缩
    img_dct_comp = dct_compress(test_image, ratio);
    mse_dct = mean((test_image(:) - img_dct_comp(:)).^2);
    psnr_dct = 10 * log10(1 / mse_dct);
    
    % SVD优势
    psnr_advantage = psnr_svd - psnr_dct;
    
    fprintf('%6.1f    %8.2f    %8.2f    %+7.2f\n', ...
            ratio, psnr_svd, psnr_dct, psnr_advantage);
end

fprintf('\n=== SVD图像压缩实验总结 ===\n');
fprintf('1. SVD提供了图像的最优低秩近似\n');
fprintf('2. 压缩质量与保留的奇异值数量直接相关\n');
fprintf('3. 不同类型图像的压缩效果差异很大\n');
fprintf('4. 自适应秩选择可以优化压缩效果\n');
fprintf('5. 分块处理可以提高计算效率\n');
fprintf('6. 彩色图像需要分通道处理\n');
fprintf('7. SVD在某些情况下优于传统压缩方法\n\n');

end

function ssim_val = compute_ssim(img1, img2)
% 简化的SSIM计算
mu1 = mean(img1(:));
mu2 = mean(img2(:));
sigma1_sq = var(img1(:));
sigma2_sq = var(img2(:));
sigma12 = mean((img1(:) - mu1) .* (img2(:) - mu2));

c1 = 0.01^2;
c2 = 0.03^2;

ssim_val = ((2*mu1*mu2 + c1) * (2*sigma12 + c2)) / ...
           ((mu1^2 + mu2^2 + c1) * (sigma1_sq + sigma2_sq + c2));
end